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Abstract 



Macroscopic models which distinguish the longitudinal and transverse temperatures can provide 
improved descriptions of the microscopic shock structures as revealed by molecular dynamics simu- 
lations. Additionally, we can include three relaxation times in the models, two based on Maxwell's 



'n , viscoelasticity and its Cattaneo-equation analog for heat flow, and a third thermal, based on the 

Krook-Boltzmann equation. This approach can replicate the observed lags of stress (which lags 



g : behind .he strain rate) and heat flux (which lags behind tire temperature gradient), as well as 

^P ■ the eventual equilibration of the two temperatures. For profile stability the time lags cannot be 

too large. By partitioning the longitudinal and transverse contributions of work and heat and 

K> ' including a tensor heat conductivity and bulk viscosity, all the qualitative microscopic features of 

strong simple-fluid shockwave structures can be reproduced. 
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FIG. 1: The steady flow shown here pictures cold material, moving to the right at speed Us and 
decelerated by slower hot fluid, moving to the right at speed Ug — Up, with Ug and Up chosen to 
fix the Shockwave location in space, Uwave = 0. An alternative way to generate such Shockwaves is 
shown in Figure 2. 

I. INTRODUCTION: PROPERTIES OF DENSE-FLUID SHOCKWAVES 

Stationary Shockwaves provide the simplest possible opportunity for the study of highly 
nonlinear transport in dense fluids. In the shock-centered steady-state coordinate frame, 
the nonequilibrium shock process converts an incoming steady stream of "cold" material 
into an outgoing stream of "hot" fluid. See Figure 1. Shockwave gradients can be huge, 
with strainrates in the terahertz range and correspondingly large pressure and temperature 
gradients, 10^^ atmospheres/centimeter and 10^^ kelvins/centimeter^ . Despite the wildly 
irreversible nature of such a nonequilibrium conversion, so long as the shock is stationary 
the overall internal energy change, E^ — Eq, can be expressed in terms of the equilibrium 
pressures and volumes of the incoming and outgoing streams of fluid: 

AE = Eh-Ec = {Ph + PcWc - Vh)/2 . 



In the steady-state coordinate frame centered on the shockwave (Figure 1), the incoming 
cold material, moving at the shock velocity Us-, is decelerated to Ug — Up by the shockfront, 
where Up is the "particle", or "piston", velocity. 

The Hugoniot relation for the energy change AE, just given, can be derived by eliminating 
the two velocities Ug and Up from the three conservation equations for mass, momentum, and 
energy^. An alternative shock-creation mechanism, quite practical for computer simulation, 
uses the symmetric collision of two blocks of cold fluid. For problems with a nonzero initial 
pressure confining pistons are required. In either case the two blocks approach each other 
with velocities ±tip, and generate two mirror-image Shockwaves identical in structure to 
those obtained with steady-state boundary conditions. See again Figure 1, as well as Figure 
2, for the geometries of these two methods for generating Shockwaves. 

Over the last forty years a wide variety of atomistic shockwave simulations, based on 
molecular dynamics, have been carried out^"— . These particle-based simulations established 
three interesting facts which simplify numerical treatments of Shockwaves. First fact: the 
boundary conditions enclosing the shockwave can be implemented easily because they are 
simply equilibrium states when viewed in a moving coordinate system. See again Figures 1 
and 2. Second fact: shockwave thicknesses are indeed only a few mean free paths^^i^A, as 
predicted for gases by numerical solutions of the Boltzmann equationi2,"iI. The small scale 
of Shockwaves makes molecular dynamics simulations relatively simple to carry out. Third 
fact: one-dimensional Shockwaves are stable^, as shown in Figure 3. Stability means that it 
is sensible to measure and compute shockwave profiles in which density, velocity, and energy 
are all expressed as functions of a single longitudinal coordinate, here chosen to define the 

CC clXlS. 

In addition to these simplifying facts there are three more facts which complicate rather 
than simplify numerical treatments. They deserve more discussion and form the heart of 
the present work: fourth fact: temperature within the shockwave is a tensor^ with different 
longitudinal and transverse values. Mott-Smith predicted the details of this complication for 
gases^-, by using an approximate bimodal velocity distribution (a spatially-varying linear 
combination of the cold and hot Maxwellian distributions). We discuss the meaning of 
"temperature" in the following Section DiL-^nlS^M.)!^, 

A fifth fact, discovered in the course of comparisons of atomistic simulations with con- 
tinuum predictions, is that the nonlocality of atomistic interactions introduces an essential 
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FIG. 2: Here two identical blocks of zero-pressure material at zizUp have collided with sufficient 
velocity to compress the fluid to twice the initial density. The two Shockwaves at the interface 
between the moving cold material and the stationary hot fluid are separating at velocities ziz{us—Up). 
The forces between particles, here and in Figures 1, 3, and 4, are short-ranged repulsive forces 
derived from the pair potential (p = (10/7r)(l — r)^. 



dependence of spatial averages on the averaging algorithm itself. Any continuum treatment 
which aims to describe two- or three-dimensional phenomena must come to grips with an ap- 
propriate choice of averaging algorithm. Lucy's one-, two-, and three-dimensional weighting 
functions used in smooth particle applied mechanicsi^"— provide a particularly appealing 



solution to the problem. Averaging is addressed in Section III. 

Last, a sixth fact, discovered more recently, is that relaxation and lag are characteristic 
of Shockwaves. Strong Shockwaves display cause-and-effect relaxation, with the shear stress, 
a = {Pyy — Prcx)/2, responding to the strain rate e = {du/dx) and the heat flux Q^ responding 
to the temperature gradient VT only after noticeable delays. These observed delay times are 
of the order of the particle-particle collision time^. Lag, relaxation, and delay are addressed 
in Section IV. 

Existing models for shockwave structure, such as the linear-transport Navier-Stokes 
equationsiii^ or the nonlinear-transport Burnett equationa^^ii^"— i2S, need to be improved 
to take these recent shock-structure observations into account. Delay has to be included 
in the models and temperature needs to have its longitudinal and transverse components 
treated separately. The present work is devoted to developing and exploring a comprehen- 
sive description of shock dynamics and developing the numerical techniques necessary to 
implement the new findings into continuum simulations. 

Following these discussions of thermal anisotropy, spatial nonlocality, and relaxation, 
we introduce a well-posed continuum model incorporating all these ideas and illustrate a 
numerical method for solving particular special cases in Section V. Section VI contains a 
summary of our results and an assessment of the prospects for future progress. 

II. KINETIC TEMPERATURE AND ITS MEASUREMENT 

Gibbs and Boltzmann related microscopic mechanics to macroscopic thermodynamics 
by showing that an ideal-gas thermometer— satisfied the Zeroth law of thermodynamic 



:sM. 



Two systems at thermal equilibrium with a Maxwell-Boltzmann ideal gas at the kinetic 
temperature 

Tgas = Teg = {pl/mk) = {pl/mk) = {pl/mk) , 

are necessarily in thermal equilibrium with each other. Thus the ideal gas is a reliable 
thermometer and can be used to measure temperature in other gases, or in liquids, or in 
solids. Consider applying the ideal-gas definition of temperature to a steady, but nonequi- 
librium, shockwave. Then there are substantial coordinate-dependent disparities between 
the longitudinal and transverse kinetic temperatures, 

{pl/mk) = T^x ; {pl/mk) = Tyy . 
5 
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FIG. 3: Six snapshots, equally spaced in time, showing the underdamped oscillation of a sinusoidal 
Shockwave. The hot shocked fluid is at twice the density of the cold unshocked material. 

These kinetic temperatures are velocity fluctuations about the local mean velocity so that 
in the comoving measurement frame the mean values of the momenta vanish: 

u{x) = (x) ; p^ = m{x - (x)) = 'm{x - u{x)) ; {p^) = {py) = . 

The kinetic definitions for the nonequilibrium longitudinal and transverse tempera- 
tures arise naturally if one imagines "measuring" them, for particular degrees of free- 
dom, by putting the nonequilibrium fluid into diagnostic contact with a comoving ideal-gas 
thermometer—. Such a thermometer is best thought of as a tiny sample of equilibrated gas, 
with the gas made up of very many very small hard particles. These thermometric particles 
undergo impulsive collisions with selected system degrees of freedom. If the ideal-gas parti- 
cles are very small the temperature measurement doesn't change the dynamical state of the 
nonequilibrium fluid^^. The ideal-gas nature of the thermometer makes it possible to ana- 
lyze the collisions from the two-body standpoint of the Boltzmann equation. Hard-disk or 
hard-sphere interactions between the thermometer and the system change, on average, the 
total kinetic energy of a system particle if it deviates from the thermometer's temperature. 
If the thermometer particles are instead pictured as parallel hard cubes (parallel squares in 
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FIG. 4: Dependence of the strainrate and the pressure tensor components Pxx and Pyy on the 
range of Lucy's weighting function w{r < h) for the stationary shockwave shown in Figure 1 (but 
with a system width eight times larger than that of Figure 1). The widths of the curves in the 
Figure increase with increasing h. 

two dimensions) with their orientations constrained, then the temperatures T^x and Tyy can 
be independently distinguished. 

The equihbrium velocity distributions, in the thermometer, are Maxwell-Boltzmann dis- 
tributions. Kinetic theory shows^ that such an ideal-gas thermometer transfers energy 
to/from a degree of freedom if the kinetic energy of that degree of freedom is less/greater 
than kTgas- When this simple mechanical definition of temperature is used to analyze shock- 
wave structure cause-and-effect relaxation and thermal anisotropy are revealed. Both these 
novel features need to be tackled and described by any realistic and comprehensive shock- 
wave model. 

III. LOCAL AVERAGES AND THEIR MEASUREMENT 



The temperature measurements just discussed require choosing a velocity for the ther- 
mometer. It must be comoving with the material in order to measure fluctuations. But 
exactly what is the velocity about which the fluctuations are measured? A useful answer 
can be found based on the weight functions used in smooth particle applied mechanics, 
"SPAM"— 122. Lucy suggested that averages, at a fixed point in space, be computed using a 



weight function w{r < h) centered there, with an arbitrary width h, and normahzed, so that 
the integral of w{r) over all space is unity. The simplest weight function is a polynomial 
chosen so that it has a smooth maximum at the origin, r = 0, and two continuous deriva- 
tives everywhere. These requirements guarantee that averages computed with the weight 
function, 

(F(r)) = Y.^rjFj/Y,Wrj ; Wrj = w{\r - rj\) , 
J j 

where F{r) is a "field variable" like density, velocity, temperature, or stress, have also two 
continuous spatial derivatives. In the sums over nearby particles {j} it is usual to choose 
the range h so that several dozen particles are included. 

These conditions on the weight function are sufficient to determine its functional form: 

WLucyir <h)oc[l- Q{r/hf + %{r/hf - 3(r//i)^] . 

Hardy's approach^^ to defining averages in Shockwaves uses the same idea as Lucy's. Evi- 
dently h must be large enough to avoid wiggles in the resulting averages, while remaining 
sufficiently small for averages to be local and inexpensive to compute. In Shockwaves a 
value for h of about three times the interparticle spacing is a good choice. Figure 4 shows 
explicitly the dependence of the pressure tensor and the velocity gradient averages on the 
range of the weight function. 

When constructing continuum models designed to reproduce atomistic simulations it is 
essential to specify the spatial averaging technique. The fact that the resulting constitu- 
tive equation depends on h is simply a reminder that atomistic mechanics and continuum 
mechanics, though similar, are not the same. 

IV. EMPIRICAL RELAXATION MODELS FOR STRESS AND HEAT FLUX 

A. Maxwell's Stress Relaxation Model and its Extension to Heat Flux 

Maxwell modeled the stress relaxation characteristic of viscoelastic fluids by introducing 
a stress relaxation time t^: 

a + T^a = r]e . 

In the Shockwave problem r] is the shear viscosity and e is the strainrate, (du/dx). Both 
time derivatives, indicated by the superior dots, are comoving with the fluid. In the ab- 
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sence of relaxation, To- = 0, Maxwell's fluid model reduces to the usual Newtonian viscous 
incompressible fluid, with shear stress a proportional to the instantaneous value of e. In the 
absence of any imposed strainrate (e = 0) the initial stress decays with a characteristic relax- 
ation time To-. For a delta-function strain rate, at t = 0, the stress has decays exponentially 
from its initial value: 

e = 5(t = 0) -^ fj = (r//r)e-*/" . 

For a relatively-simple case, with t] = t = 1, and a localized strain rate, like that in the 
Landau-Lifshitz description of a weak shock-: 

1 



Maxwell's model has an analytic solution: 



a 



(t) = e~* In Vl + e+2* . 



Figure 5 illustrates the stress response for the Newtonian case r = 0, and for two Maxwellian 
relaxation times, t„ = 1 and t„ = 4. Recent molecular dynamics shockwave simulations have 
shown that both stress and heat flux exhibit delayed responses^. 

Exactly the same ideas can be, and have been^, applied to heat flux. If we introduce the 
relaxation time tq into Fourier's law for heat flow, the result is the Cattaneo equation: 

and the heat flux lags behind the temperature gradient by a time of the order of r. 

On physical grounds the time derivative here is again comoving with the fluid. The Cat- 
taneo equation describes heat flux and predicts its decay, just as did Maxwell's formulation 

of stress decay: 

1 



VT ot _, , ^+f — > Qif) oc -e"* In Vl + e+2* . 



In the following Section we illustrate how to incorporate these relaxation effects for stress 
and heat flux into a simple dense-fluid shockwave model. 

B. Krook-Boltzmann Thermal Relaxation 

The Boltzmann equationi^ models the dynamics of a dilute gas in which the gas particles 
undergo occasional two-body collisions. The evolution of the velocity distribution function. 



0.5 - 
0.4 - 
0.3 - 
0.2 - 
0.1 - 
O.O 



Maxwell Stress Relaxation 

T = O 




-10 < t < +10 

FIG. 5: Delayed stress a in response to the strain rate l/[e~* + e~^^] with unit viscosity, rj = 1. 
Maxwell's relaxation time r controls the stress response: a + t& = ije. 

f{p,r,t) for the phase-space density of particles with momentum p at location r at time t, 
can be approximated by an exponential relaxation toward equilibrium: 

(df/dt) ^ [/e, - f]/T ^^ / + Tidf/dt) = /e, . 

This approximate "Krook-Boltzmann" equation has exactly the same form as do the Maxwell 
and Cattaneo relaxation equations. Here the relaxation time r defined by this approximate 
equation is of the order of the mean collision time. Because two-body "conservative" colli- 
sions conserve mass, momentum, and energy, the equilibrium distribution, toward which / 
relaxes, necessarily has the same density, stream velocity, and energy, as does the nonequi- 
librium distribution /. 

Shockwaves convert macroscopic longitudinal kinetic energy into microscopic "thermal" 
internal energy, 

AmV2 — > Ae , 

through collisions, so that it is reasonable to expect shockwave stresses and temperatures 
to relax and equilibrate in a time of order the collision time r. We include these delay and 
relaxation effects in the macroscopic model formulated in the next Section. 
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V. FORMULATION OF A MACROSCOPIC MODEL 

Any solution of the continuum evolution equations, 

p = -pV -u; pu = -V ■ P ; pe = -Vm : P - V ■ Q , 

requires constitutive models giving the pressure tensor P and heat flux Q in terms of the 
underlying variables { p,u,e }, the density, velocity, and energy per unit mass. For com- 
pleteness, in view of the relaxational results from molecular dynamics simulations, we must 
include separate tensor temperature components, T^j.^ and Tyy, in the list of state variables. 
In a gas the difference is simply related to the pressure tensor: 

y-'- XX yy) \^xx -^yy ) / \P'^ ) 5 

where k is Boltzmann's constant per unit mass. 

In a dense fluid, the potential contribution to anisotropicity is comparable to the kinetic 
contribution^^. A semiquantitative description of the potential part of the shear stress results 
if the equilibrium fluid structure is sheared, at the strainrate e, for the Maxwell relaxation 
time To-. The shear distortion of the pair distribution function in dense fluids has been 
studied experimentally^ and modeled with molecular dynamics^. In both cases Maxwell's 
relaxation provides a good description of the potential contribution to the shear stress. Thus 
the gas-phase description of shear anisotropy must be modified in order to describe dense 
fluids. 

To solve this problem we choose to separate the work and heat contributing to energy 
change into separate longitudinal and transverse parts. The simplest choice is a time- 
independent division of work and heat into longitudinal and transverse parts: 

-aVu : P -^ AT,, ; -(1 - a)Vu : P ^ ATyy ; 

-/3V • Q ^ AT,, ; -(1 - /3)V ■ g ^ ATyy . 

The Navier-Stokes equations correspond to the choice a = f3 = 1/2. In a shockwave, where 
the kinetic energy is initially longitudinal, we would expect instead a ^ /3 ~ 1. 

To explore the consequences of this division we consider in what follows a simple van 
der Waals model, with the energy and equilibrium pressure expressed as sums of density- 
dependent and temperature-dependent contributions. A slightly more flexible model^ can 
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be based on Griineisen's separation of the energy and pressure into corresponding "cold" 
and "thermal" parts. 

Away from equilibrium we include Maxwell's delayed viscous response in the pressure 
tensor. For a two-dimensional fluid undergoing uniaxial compression and with shear viscosity 
rj and vanishing bulk viscosity, we have: 

Pxx = Peq - cr ; Pyy = Pcci + <J ; a + T^a = r]{du/dx) . 

The stress relaxation time Tq- describes the delay in the response of the shear stress a to the 
strainrate e = (du/dx). 

If the longitudinal and transverse temperatures are constrained to differ, we would expect 
the stationary nonequilibrium heat flux vector to obey a tensor form of Fourier's law: 

Kq^X '^XX\^'^ XX/ ^•^ ) riyy\JXl yy I (XX j . 

In the Shockwave problem the effects of delay and eventual equilibration both need to be 
included. For simplicity we add on corresponding delays and thermal relaxation to the 
continuum evolution equations for the heat flux and the temperatures: 

Qx 3 ~QxnQ ] Txx 3 [Tyy — T^x) / ^T ] Tyy D [Txx " Tyy) / Tj^ ■ 

To model a dense A^-particle van der Waals fluid, as opposed to a dilute gas, we approx- 
imate the potential part of the thermal energy by setting it equal to the kinetic part: 

-E'* — -Ecold — Ek = Nk{Txx + Tyy)/ 2 )■ -^Thermal = Nk{Txx + Tyy) . 

The motivation for studying such simple continuum models derives from the results of 
molecular dynamics simulations of stationary Shockwaves^"—. Just as in the continuum case, 
these microscopic molecular dynamics simulations conserve mass, momentum, and energy, 
so that the stationary fluxes of these quantities, 

pu , Pxx + pu^ ; pu[e + [PxxIp) + (mV2)] + Qx , 

are constant throughout the flow. These simulations show further that both the stress 
and heat flux lag behind the strainrate and temperature gradient. The lags are physically 
reasonable from the collisional cause-and-effect standpoint. 
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Newton's viscosity and Fourier's heat conduction both describe instantaneous relation- 
ships. Taken hterally, these two hnear laws imply that the stress a and heat flux Q respond 
instantaneously, and supersonically, to the strainrate e and temperature gradient VT. 

Certainly such an instantaneous response is impossible. If we imagine reversing a time- 
reversible Newtonian motion of the shock process, another apparent shortcoming of the 
Navier-Stokes formulation is revealed. Newton's and Fourier's laws, 

cr oc e ; Qx oc —{dT/dx) , 

if applied to a time-reversible flow, imply that the stress changes sign (as e changes sign 
when the motion is reversed) while the heat flux does not (as the temperature gradient 
has no time-dependence). Both conclusions are inconsistent with time- reversible Newtonian 
dynamics. 

A detailed atomistic analysis of the pressure tensor and the heat flux vector^^ shows 
that these functions are respectively even and odd functions of time, so that Newton's and 
Fourier's ideas are necessarily inexact as they lack the proper delay time inherent in inter- 
particle collisions. Lacking a more fundamental approach to time-reversible irreversibility, 
we seek to learn more by exploring explicitly the irreversible nature of continuum models. 

VI. NUMERICAL SOLUTIONS OF THE VAN DER WAALS SHOCKWAVE 
MODEL 

A one-dimensional "staggered grid", with density evaluated within Nc cells of length 
dx, bounded by Nn = Nc + 1 nodes, and with the remaining long list of time-dependent 
variables {u^e^T^xiTyy^a^Q} given at the nodes, provides a basis for an iterative solution 
of the continuum equationsi^i2^. Our assumption relating the energy change to the changes 
in the two temperatures gives the set of nodal variables {m, T^x-, Tyy, cr, Q} with the internal 
energy density given by 

e = (p/2) + kTxx + kTyy . 

We have found that such an approach can be applied to wave-structure relaxation with 
longitudinal and transverse work and heat separation, as well. The Landau-Lifshitz weak- 
shock solution^ - for constant shear viscosity and thermal conductivity, and without any 
relaxation - makes a useful initial condition. In the stationary-shockwave coordinate system 
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FIG. 6: Development of the stationary temperature profiles in a Shockwave with shear viscosity 

{rj = 4), heat conductivities {k^x = Hyy = 2} and relaxation times {t^ = tq = Tr = 1} ■ The 
converged temperatures are shown at the right. Here the work done and the heat transfer initially 
affect only the longitudinal temperature. 

errors in the initial condition move to the right, away from the shockwave, at approximately 
the speed of sound. See Figures 6 and 7 for transient results from typical solutions of the 
continuum equations.. 

A successful numerical evolution algorithm proceeds from an initial guess by iterating a 
series of four steps: (i) specify the six dependent variables {pc, Vn,Txxn,Tyyn, cTn, Qn} at all 
the interior cells and nodes; (ii) compute all the remaining variables and all the gradients 
with centered sums and differences: 

Uc{x) = [un{x — dx/2) + Un{x + dx/2)]/2 ; Pn{x) = [pc{x — dx/2) + pc{x + dx/2)]/2 ; 

dTii/dx = [Tii{x + dx/2) — Ta^x — dx/2)]/dx ; 

(iii) compute the righthandsides of the five sets of differential equations. 

For instance, the change in energy at a particular node could be evaluated as follows: 

{de/dt)n = -Un{de/dx)n - [{Pxxdu/dx)n + {dQx/dx)n]/pn ; 

(iv) use the fourth-order Runge-Kutta method to integrate the A''c + 5A^„ ordinary differential 
equations for one timestep dt, providing the information necessary for a return to step (i) 
for the execution of the next timestep. The numerical values of the mass, momentum, and 
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FIG. 7: Development of the stationary temperature profiles in a shockwave with shear viscosity 

{rj = 4), heat conductivities {kxx = Hyy = 2} and relaxation times {t^ = tq = r^ = 1} . The 
converged temperatures are shown at the right. Here the work done and the heat transfer initially 
affect only the transverse temperature. 

energy, as well as their fluxes can be used to help estimate the initial conditions. For the 
twofold compression shockwave we use to illustrate these ideas, the fluxes and boundary 
values are the following: 



eq 



pe ; e = (p/2) + T^x + Tyy ; Ceq = (p/2) + 2T ; 



pu = 2; Pxx + pu^ = 9/2 ; {pu)[e + {Pxxl p) + (m'/S)] + Qa= = 6 ; 
p : (2 ^ 1) ; M : (1 ^ 2) ; P : (1/2 ^ 5/2) ; T : (0 ^ 1/8) ; e : (1/2 ^ 5/4) . 

Here the cold and hot boundary values are linked by arrows: {cold — )■ hot). Both Qx and a 
necessarily vanish at the boundaries, Q^ : (0 — )■ 0) ; a : (0 — )> 0) . 

Figures 6, 7, and 8 illustrate typical solutions. In order to circumvent numerical in- 
stabilities in the numerical work one can (i) increase the number of cells, (ii) reduce the 
timestep and/or cell size, (iii) introduce an explicit artificial time- dependence in the param- 
eters {r], K, r} in order to enhance convergence. In this way we have obtained solutions of 
the continuum shockwave model for a wide range of parameters. The same ideas can be used 
to study special cases in which stress or heat flux are not delayed or in which temperature is 
scalar rather than tensor. The sample solutions shown in Figures 6-8 show how the partition 
of heat and work can affect the stationary shockwave. 
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FIG. 8: Stationary profiles showing the longitudinal and transverse temperature gradients 
[{dTxx/dx) is solid; [dTyy/dx) is dashed] as well as the velocity gradients. Notice that the heat flux 
and shear stress lag behind the gradients which "cause" them. 

VII. CONCLUSIONS AND PROSPECTS 

The simulation of nonequilibrium stationary states with molecular dynamics, particularly 
in the two-block geometry of Figure 2, emphasizes Loschmidt's reversibility paradox— i^i^i^. 
Evidently a movie of exactly the same dynamical states, played backward in time, satisfies 
all the microscopic motion equations. Such reversed motions contradict macroscopic physics 
and are never observed in practice. They would be inherently Lyapunov unstable, with any 
small perturbation (such as roundoff in the last place) growing exponentially in time and so 
destroying the reversed trajectory. 

Because time-reversed solutions of the Newtonian equations of motion are not observable 
it is legitimate to use time-irreversible models in interpreting the solutions. The noticeable 
time-delays, for both stress and heat flux, observed in these solutions legitimates also the use 
of Maxwell-Cattaneo-Krook relaxation. These innovations are useful to the goal of flnding 
macroscopic descriptions conforming to microscopic observations. 

Although we have been able to flnd stable solutions for the most general description 
considered here (three relaxation times, partition of heat and work, tensor temperature) 
there are stringent limits on the parameter ranges for which such solutions exist. On physical 
grounds stress and heat flux relaxation must be relatively rapid. A fluid's memory cannot be 
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too long. A systematic study of stability is complicated by the large number of parameters 
involved. Nevertheless, carefully chosen example cases should shed additional light on the 
physics of relaxation and of strong Shockwaves. At the moment the step of generalizing the 
physical ideas further, for instance by considering the state dependence of the relaxation 
times, is premature. But we can confidently expect progress there in the future. 
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